library(rms); library(margins)
load("data-2006-cces.RData")

# Subset to only respondents that correctly identified the PID of the Senator they were evaluating
only.pid.correct <- data[which(data$sen_pid_correct == 1),]

###########
# Table 1 #
###########
model.1 <- glm(sen_approval_binary ~ state_pid + upforelection + balance + opp_party + ind + decades + factorpid + female 
               + black + latino + other + college + age, data=only.pid.correct, family=binomial(link = "logit"))

model.2 <- glm(sen_approval_binary ~ state_pid + upforelection + balance*opp_party + opp_party + balance*ind + ind + decades 
               + factorpid + female + black + latino + other + college + age, data=only.pid.correct, 
               family=binomial(link = "logit"))

# Calculating clustered SEs (by respondent) for Table 1 models
cluster.model.1 <- lrm(sen_approval_binary ~ state_pid + upforelection + balance + opp_party + ind + decades + factorpid 
                       + female + black + latino + other + college + age, data=only.pid.correct, x=T, y=T)
robcov(cluster.model.1, cluster=only.pid.correct$id)

cluster.model.2 <- lrm(sen_approval_binary ~ state_pid + upforelection + balance*opp_party + opp_party + ind*balance + decades 
                   + factorpid + female + black + latino + other + college + age, data=only.pid.correct, x=T, y=T)
robcov(cluster.model.2, cluster=only.pid.correct$id)

############
# Figure 2 #
############
# Colors for figure
mygray <- rgb(col2rgb("gray80")[1], col2rgb("gray80")[2], col2rgb("gray80")[3], alpha = 100, maxColorValue = 255) 

# Generating predicted values, holding covariates at their mean values
# First 13 predicted values are for independents (ind=1). Second 13 are for copartisans
independent.dataset <- data.frame(balance=rep(c(-.5,-.4,-.3,-.2,-.1,0,.1,.2,.3,.4,.5,.6,.7),2),
                                  state_pid=rep(mean(only.pid.correct$state_pid), 26),
                                  upforelection=rep(mean(only.pid.correct$upforelection),26),
                                  opp_party=c(rep(0,13),rep(0,13)),
                                  ind=c(rep(1,13),rep(0,13)),
                                  decades=rep(mean(only.pid.correct$decades),26),
                                  factorpid=as.factor(rep(2,26)),
                                  female=rep(mean(only.pid.correct$female),26),
                                  black=rep(mean(only.pid.correct$black),26),
                                  latino=rep(mean(only.pid.correct$latino),26),
                                  other=rep(mean(only.pid.correct$other),26),
                                  college=rep(mean(only.pid.correct$college),26),
                                  age=rep(mean(only.pid.correct$age),26))
# First 13 predicted values are for copartisans. Second 13 are for opposite partisans (opp_party=1).
opp.party.dataset <- data.frame(balance=c(-.5,-.4,-.3,-.2,-.1,0,.1,.2,.3,.4,.5,.6,.7),
                                state_pid=rep(mean(only.pid.correct$state_pid), 26),
                                upforelection=rep(mean(only.pid.correct$upforelection),26),
                                opp_party=c(rep(0,13),rep(1,13)),
                                ind=rep(0,26),
                                decades=rep(mean(only.pid.correct$decades),26),
                                factorpid=as.factor(rep(2,26)),
                                female=rep(mean(only.pid.correct$female),26),
                                black=rep(mean(only.pid.correct$black),26),
                                latino=rep(mean(only.pid.correct$latino),26),
                                other=rep(mean(only.pid.correct$other),26),
                                college=rep(mean(only.pid.correct$college),26),
                                age=rep(mean(only.pid.correct$age),26))

opp.party.fitted <- predict(model.2, opp.party.dataset, type=c("response"), se.fit=T)
independent.fitted <- predict(model.2, independent.dataset, type=c("response"), se.fit=T)

opp.marginals <- opp.party.fitted$fit[1:13] - opp.party.fitted$fit[14:26]
opp.ses <- sqrt((opp.party.fitted$se[1:13])^2 + (opp.party.fitted$se[14:26])^2)

independent.marginals <- independent.fitted$fit[14:26] - independent.fitted$fit[1:13]
independent.ses <- sqrt((independent.fitted$se[1:13])^2 + (independent.fitted$se[14:26])^2)

opp_low_ci<- opp.marginals-1.96*opp.ses
opp_high_ci<-opp.marginals+1.96*opp.ses
balance_values<-c(seq(-.5, 0.7, by=.1))

ind_low_ci<- independent.marginals-1.96*independent.ses
ind_high_ci<- independent.marginals+1.96*independent.ses
balance_values<-c(seq(-.5, 0.7, by=.1))

# Making Figure 2
par(mar=c(3, 2.25, 1, 1), oma=c(0,.5,0,0))
plot(balance_values, opp.marginals, type="n",
     xlim=c(-.5,.7), ylim=c(0,.8),
     xaxt="n", yaxt="n", xlab="", ylab="",
     main="",
     cex=0.9)
polygon(c(balance_values,rev(balance_values)),
        c(opp_low_ci,rev(opp_high_ci)), col=mygray, border=NA)
lines(balance_values, opp.marginals)

polygon(c(balance_values,rev(balance_values)),
        c(ind_low_ci,rev(ind_high_ci)), col=mygray, border=NA)
lines(balance_values, independent.marginals, lty=2)

rug(data$balance, ticksize = 0.03, side = 1, lwd = 0.5, col = par("fg"),
    quiet = getOption("warn") < 0)

axis(1, tck=0, cex.axis=0.8, cex.lab=0.8, mgp=c(0.5, 0.1, 0), lty=0, labels=TRUE)
axis(2, tck = 0,las=2, lwd = 0, cex.axis = 0.8, cex.lab = 1, mgp=c(0.5, 0.1, 0), 
     lty=0)
title(xlab = "Policy emphasis                                            Non-policy emphasis", line = 1.25, cex.lab = 1.15, ylab="Difference in predicted probability of Senator approval")
text(.2, .7, "Copartisans vs. opposite partisans",cex=1)
text(-.1, .3, "Copartisans vs. Independents",cex=1)

############
# Figure 3 #
############
model.f3 <- glm(sen_approval_binary ~ state_pid + upforelection + balance*opp_party*pid_align + opp_party + decades 
                + factorpid + female + black + latino + other + college + age, data=only.pid.correct, 
                family=binomial(link = "logit"))
summary(model.f3)

# Generating predicted values, holding covariates at their mean values
pid.not.aligned.dataset <- data.frame(balance=rep(c(-.5,-.4,-.3,-.2,-.1,0,.1,.2,.3,.4,.5,.6,.7),2),
                                      state_pid=rep(mean(only.pid.correct$state_pid), 26),
                                      upforelection=rep(mean(only.pid.correct$upforelection),26),
                                      opp_party=c(rep(0,13),rep(1,13)),
                                      decades=rep(mean(only.pid.correct$decades),26),
                                      factorpid=as.factor(rep(2,26)),
                                      female=rep(mean(only.pid.correct$female),26),
                                      black=rep(mean(only.pid.correct$black),26),
                                      latino=rep(mean(only.pid.correct$latino),26),
                                      other=rep(mean(only.pid.correct$other),26),
                                      college=rep(mean(only.pid.correct$college),26),
                                      age=rep(mean(only.pid.correct$age),26),
                                      pid_align=rep(0,26))
pid.aligned.dataset <- data.frame(balance=c(-.5,-.4,-.3,-.2,-.1,0,.1,.2,.3,.4,.5,.6,.7),
                                  state_pid=rep(mean(only.pid.correct$state_pid), 26),
                                  upforelection=rep(mean(only.pid.correct$upforelection),26),
                                  opp_party=c(rep(0,13),rep(1,13)),
                                  decades=rep(mean(only.pid.correct$decades),26),
                                  factorpid=as.factor(rep(2,26)),
                                  female=rep(mean(only.pid.correct$female),26),
                                  black=rep(mean(only.pid.correct$black),26),
                                  latino=rep(mean(only.pid.correct$latino),26),
                                  other=rep(mean(only.pid.correct$other),26),
                                  college=rep(mean(only.pid.correct$college),26),
                                  age=rep(mean(only.pid.correct$age),26),
                                  pid_align=rep(1,26))

aligned.fitted <- predict(model.f3, pid.aligned.dataset, type=c("response"), se.fit=T)
not.aligned.fitted <- predict(model.f3, pid.not.aligned.dataset, type=c("response"), se.fit=T)

aligned.marginals <- aligned.fitted$fit[1:13] - aligned.fitted$fit[14:26]
aligned.ses <- sqrt((aligned.fitted$se[1:13])^2 + (aligned.fitted$se[14:26])^2)

not.aligned.marginals <- not.aligned.fitted$fit[1:13] - not.aligned.fitted$fit[14:26]
not.aligned.ses <- sqrt((not.aligned.fitted$se[1:13])^2 + (not.aligned.fitted$se[14:26])^2)

aligned_low_ci<- aligned.marginals-1.96*aligned.ses
aligned_high_ci<-aligned.marginals+1.96*aligned.ses
balance_values<-c(seq(-.5, 0.7, by=.1))

not_aligned_low_ci<- not.aligned.marginals-1.96*not.aligned.ses
not_aligned_high_ci<- not.aligned.marginals+1.96*not.aligned.ses

# Making Figure 3
par(mar=c(3, 2.25, 1, 1), oma=c(0,.5,0,0))
plot(balance_values, aligned.marginals, type="n",
     xlim=c(-.5,.7), ylim=c(0,.9),
     xaxt="n", yaxt="n", xlab="", ylab="",
     main="",
     cex=0.9)
polygon(c(balance_values,rev(balance_values)),
        c(aligned_low_ci,rev(aligned_high_ci)), col=mygray, border=NA)
lines(balance_values, aligned.marginals, lty=2)

polygon(c(balance_values,rev(balance_values)),
        c(not_aligned_low_ci,rev(not_aligned_high_ci)), col=mygray, border=NA)
lines(balance_values, not.aligned.marginals)

rug(data$balance, ticksize = 0.03, side = 1, lwd = 0.5, col = par("fg"),
    quiet = getOption("warn") < 0)

axis(1, tck=0, cex.axis=0.8, cex.lab=0.8, mgp=c(0.5, 0.1, 0), lty=0, labels=TRUE)
axis(2, tck = 0,las=2, lwd = 0, cex.axis = 0.8, cex.lab = 1, mgp=c(0.5, 0.1, 0), 
     lty=0, at=c(0.1,0.3,0.5,0.7,0.9))
title(xlab = "Policy emphasis                                            Non-policy emphasis", line = 1.25, cex.lab = 1.15, ylab="Difference in predicted probability of Senator approval")
text(.2, .78, "Ideology aligned with partisanship",cex=1)
text(-.15, .45, "Ideology not aligned with partisanship",cex=1)
